home *** CD-ROM | disk | FTP | other *** search
/ SGI Hot Mix 17 / Hot Mix 17.iso / HM17_SGI / research / lib / obsolete / kmeans.pro < prev    next >
Text File  |  1997-07-08  |  13KB  |  491 lines
  1. ; $Id: kmeans.pro,v 1.2 1997/01/15 04:02:19 ali Exp $
  2. ;
  3. ;  Copyright (c) 1991-1997, Research Systems Inc.  All rights
  4. ;  reserved. Unauthorized reproduction prohibited.
  5.  
  6.  
  7.  
  8. Function Normal1, Data,R,C
  9.  
  10. ;Normal1 returns the matrix obtained by normalizing the columns
  11. ; of Data.  R= # of rows, C= # of columns.
  12.  
  13. Y= Data-Data#Replicate(1./R,R) #Replicate(1.,R)
  14. std =sqrt(Y^2 # Replicate(1./(R-1),R))
  15. D1=Fltarr(c,c)
  16. for i=0l,C-1 do              $
  17.  if std(i) NE 0 then D1(i,i)=1./std(i) else d1(i,i)=0
  18.  
  19. return, D1#Y
  20. end     ;Normal1
  21.  
  22.  
  23.  
  24.  
  25.  
  26.  Pro StatComp,D,V,N,Mx,Mn,STD,SS
  27.  
  28.  ; compute summary statistics,min,max,standard deviation for
  29.  ; all variables in a cluster whose cases are listed in V. 
  30.  ; Statistics stored  in the 1-dim  arrays, Mx=Max,Mn=Min,
  31.  ; STD=standard deviation,SS=sum  of squares. D=Data,  see
  32.  ; procedure kmeans.
  33.  
  34.  R=Size(D)               
  35.  N=R(1)                    ; N = # of variables
  36.  R1=R(0)               
  37.  R=R(2)                    ; R = # of cases
  38.  
  39.  
  40.  Y= D-V#Replicate(1.,R)    ; Compute sum of squares and
  41.                            ; standard deviation
  42.  if(R1 GT 1) THEN SS= Y^2#Replicate(1.,R) ELSE SS=Y^2
  43.  
  44.  if(R1 GT 1) THEN  STD=SQRT(Y^2#Replicate(1./(R-1),R))     $
  45.   ELSE STD=FLtArr(N)
  46.  
  47.  Mx=Fltarr(N)
  48.  MN=Mx
  49.  
  50.  for i=0l,N-1 DO BEGIN      ; compute max and min
  51.      MN(i)=min(D(i,*))
  52.      Mx(i)=max(D(i,*))
  53.  ENDFOR
  54.  RETURN
  55.  END
  56.   
  57.  
  58.  
  59.  
  60.  
  61.  
  62. Pro VAnova,Data,Cluster,VarNames,SS,unit,N,CN,R,SizeB,B1,SW1
  63.  
  64. ;Compute and print out 1-way  analysis of variance for each
  65. ;variable. Interpret clusters as treatments. N=# of variables,
  66. ; CN= #of clusters, R=of cases, SizeB(i)= size of ith cluster.
  67.  
  68.  
  69. Common KBlock,B                  ;B(i,j)= ith variable mean
  70.                                  ; for cluster j. 
  71.  
  72. printf,unit,Format='(/,/,/)'
  73. printf,unit,"                    Overall Variable Statistics"
  74. printf,unit,"*******************************************************************"
  75.  
  76.   printf,unit,                 $
  77. Format='("Variable",6x,"Within SS",4X,"DF",4x,"Between SS",4X,"DF",5X,"FRatio")'
  78.  
  79. SW1=Fltarr(N)
  80. MS= Replicate(1.,N)#SizeB       ;compute between sum of squares
  81. TT= B*MS#Replicate(1.,CN)       
  82. TT= TT*TT/Total(SizeB)
  83. B1= B^2 * MS#Replicate(1.,CN)-TT  ; B1(i) = between SS for ith 
  84.                                   ;variable 
  85. DFB=CN-1                          ; compute between degrees of freedom 
  86. DFW=R-DFB-1                  ; compute within degrees of freedom
  87. if DFB EQ 0 THEN DFB= 1.e-10
  88. if DFW EQ 0 THEN DFW = 1.e-10
  89.  
  90. for i=0l,N-1 DO BEGIN
  91.     SSB=B1(i) 
  92.     SW= Total(SS(i,*))            ; within SS
  93.     SW1(i)=SW
  94.   printf,unit,Format='(A8,4X,G11.4,3X,I3,3X,G11.4,3X,I3,1X,G11.4)',VarNames(i),SW,DFW,SSB,DFB,SSB*(DFW)/(SW*DFB)
  95. ENDFOR
  96.  
  97.  RETURN
  98. END                 ;Procedure Vanova
  99.     
  100.     
  101.  
  102.  
  103.  
  104.  
  105.  Pro OutPut1, Clus,Data,CaseNam,VarNam,unit,SD1,SB1,SW1
  106.  ; Output mean,min,max and standard deviation for each
  107.  ; variable for each cluster.Also, output overall analysis
  108.  ; of variance.
  109.  
  110.  Common KBlock,B 
  111.         ;B(i,j) = value of ith variable in jth cluster
  112.  
  113.  SB=Size(B)
  114.  K= SB(2)                ; final # of clusters
  115.  N= SB(1)                ; number of variables
  116.  SD=Size(Data)
  117.  CN=SD(2)                ; number of cases
  118.  SD1=Fltarr(N,K)
  119.  NC=Size(CaseNam)
  120.  
  121.  if(NC(1) NE 0 ) THEN BEGIN 
  122.                  ;check and, maybe, fix case names
  123.      CaseName=CaseNam
  124.      if(NC(1) LT CN) THEN BEGIN
  125.      printf,unit,'kmeans-missing Case names'
  126.      I=Indgen(CN)
  127.      CaseName=[CaseName,'Case'+STRTRIM(I(NC(1):CN-1),2.)]
  128.    ENDIF
  129.  ENDIF ELSE CaseName= 'Case'+STRTRim(INDgen(CN),2)
  130.  
  131.  NC=Size(VarNam)
  132.  
  133.  if(NC(1) NE 0 ) THEN BEGIN
  134.           ;likewise, for variable names
  135.    VarName=VarNam
  136.    if(NC(1) LT N) THEN BEGIN
  137.      printf,unit,'kmeans-missing Variable names'
  138.      I=Indgen(N)
  139.      VarName=[VarName,'Var'+STRTRIM(I(NC(1):N-1)),2]
  140.    ENDIF
  141.  ENDIF ELSE VarName='Var'+ STRTrim(INDgen(N),2)
  142.  
  143.  
  144.    
  145.  SSW=Fltarr(N,K)
  146.                         ;SSW(i)= sum of squares for cluster i
  147.  SizeB= FltArr(K)          ;SizeB(i) = # of cases in cluster i
  148.  
  149.  for i = 0l,K-1 DO BEGIN    ;print out cases and statistics for
  150.                            ;each cluster
  151.  V=B(*,i)    
  152.                  ;V= vector of mean var values for ith cluster
  153.  printf,unit,Format='(/,/,/)'
  154.  
  155.  printf,unit,Format='(20X,"Cluster:",I12)',i
  156.  
  157. printf,unit,'********************************************************************************'
  158.  
  159.  
  160.  printf,unit,                       $
  161. Format='(/,6X,"Case",5X,"Distance",3X,"|",5X,"Var",8X,"Min",8X,"Max",7X,"Mean",8X,"STD",/,/)'
  162.  
  163. CL=where(Clus EQ i,count)    
  164.                            ;Cl= list of cases in ith cluster
  165. if count ne 0 THEN BEGIN
  166.  SizeB(i)=count  
  167.  M=count> N
  168.  D1=Data(*,CL)                ;D1 is Data restricted to the
  169.                               ;cases in cluster i.
  170.  
  171.  StatComp,D1,V,N,mx, mn,STD,SS        ; compute var stats 
  172.  SSW(*,i)=SS          
  173.                       ;compute within sum of squares for anova
  174.  SD1(*,i)=STD
  175.  
  176.  for j=0l,M-1 DO BEGIN                 
  177.  
  178.      if(J LE count-1) THEN BEGIN   ;print Cases and distances
  179.        Dist= Data(*,CL(J))-V
  180.        printf,unit,                $
  181. Format='(A10,2X,G11.4,3X,A1,$)',   $
  182. CaseName(CL(j)),SQRT(Total(Dist*Dist)),"|"
  183.      ENDIF ELSE printf,unit,Format='(26X,"|",$)'
  184.  
  185.  
  186.  
  187.    if(J LE N-1) THEN BEGIN
  188.                   ; print variables and their stats
  189.  
  190. printf,unit,Format='(2X,A6,G11.4,G11.4,G11.4,G11.4)',    $
  191.                    VarName(j),mn(j),mx(j),V(j),STD(j)
  192.    ENDIF ELSE printf,unit," "
  193.  ENDFOR
  194. ENDIF
  195. ENDFor
  196.  
  197.  VAnova,Data,Clus,VarName,SSW,unit,N,K,CN,SizeB,SB1,SW1
  198.                                              ; output anova
  199.  
  200. RETURN
  201. END        ; OutPut1
  202.      
  203.  
  204.  
  205.  
  206.  
  207.   
  208. Pro KMeans1, Data,Clus,iter,K
  209.  
  210. ; Data is a (C,R) dimensioned array of R cases and C
  211. ; variables.Clus must be a one-dimensional array of size
  212. ; R and with each element between 0 and K-1. K is the number
  213. ; of clusters and iter the maximum #of times to relocate case
  214. ; in clusters to minimize error. On exit, CLUS represents
  215. ; the final clusters of the partition. Clus(i)=cluster
  216. ; containing casei.
  217.  
  218.  
  219.  
  220. Common KBlock,B 
  221.  On_Error,2  
  222.  
  223.  SD=Size(Data)
  224.  if SD(0) eq 1 then R = 1 else R=SD(2)                          ; R= # of cases
  225.  C=SD(1)                          ; N= # of variables
  226.  
  227.  B=Fltarr(C,K)                    
  228.  D=Fltarr(K,R)
  229.  INCLUS=FLtarr(R,k)
  230.  L=Fltarr(K)                      ;L(i)= size of ith cluster
  231.  
  232.  
  233.  for i=0l,k-1 DO BEGIN
  234.     X=( CLUS EQ i)
  235.     L(i)= Total(x)
  236.     if(L(i) NE 0) THEN X=X/L(i)
  237.     INCLUS(*,i)=X
  238.   END
  239.  
  240.  B= Data#INCLUS
  241.             ;B(N,L)= mean of the Nth variable in cluster L
  242.      
  243.  CD=Fltarr(R) 
  244.  
  245.  
  246.  X=Data-B(*,Clus)
  247.  e = Total(X*X)
  248.  
  249.  for j=1l,iter DO BEGIN         ; reiteratively, relocate
  250.                                ; cases in clusters
  251.                                ; to reduce error
  252.     GotOne=0
  253.     for i= 0l,R-1 DO BEGIN      ; run thru cases
  254.          Min=1.e30
  255.          NK=-1
  256.          D1= Data(*,i)-B(*,Clus(i))      ; 
  257.      if(L(CLUS(i)) GT 1)               $
  258.      THEN D1= L(CLUS(I))*Total(D1*D1)/(L(Clus(i))-1)    $
  259.      else D1=0
  260.                ;D1=~Distance from i to its cluster Clus(i) 
  261.          for n= 0l,k-1 DO BEGIN
  262.              if(n NE CLUS(i)) THEN BEGIN
  263.                 D2= Data(*,i)-B(*,n)  
  264.                                 
  265.                 D2= L(n)*Total(D2*D2)/(L(n)+1) 
  266.                        ;D2=~Distance from case i to cluster n
  267.  
  268.                 if (D2-D1 LT Min) THEN BEGIN
  269.                         ;D1-D2= error change if object i is
  270.                         ; relocated into cluster n
  271.                     Min=D2-D1
  272.                     NK=n
  273.                     ENDIF
  274.               ENDIF
  275.            ENDFOR
  276.  
  277.           if(Min LT 0) THEN BEGIN
  278.                   ;relocate to cluster NK if it reduces error
  279.         
  280.             B(*,CLus(i))=   $
  281.            (B(*,Clus(i))*L(CLus(i))-Data(*,i))/(L(Clus(i))-1)
  282.                         ;update cluster means after relocation
  283.             B(*,NK)=(B(*,NK)*L(NK)+Data(*,i))/(L(NK)+1)
  284.             L(CLUS(I))=L(CLUS(I))-1      ; update cluster size
  285.             CLUS(I)=NK
  286.             L(NK)=L(NK)+1
  287.             e=e+Min                      ; update error
  288.             GotOne=1
  289.           ENDIF
  290.   
  291.      ENDFOR
  292.  
  293.    if(Not GotOne ) THEN RETURN     ;quit with final Partition if
  294.                                  ;cant relocate to reduce error  
  295.     ENDFOR
  296.  
  297.  
  298.  
  299. RETURN
  300. END
  301.              
  302.             
  303.              
  304.      
  305.   
  306.  
  307.  
  308.  
  309. Pro kmeans,Data1,CLuster,VarName = VarName,CaseName = CaseName,$
  310.           Iter=IT,Number=n,Norm=Nr,Missing=M, List_Name=LN,$
  311.           ClustMeans=CM,ClustSTD=CS,SSBetween=SB,SSWithin=SW,NoPrint=NP
  312.  
  313.  
  314. ;+
  315. ;
  316. ; NAME:
  317. ;    KMEANS
  318. ;
  319. ; PURPOSE:
  320. ;    To split the cases in Data into n (default=2) groups that have 
  321. ;    minimum in  group variation relative to between group variation.
  322. ;
  323. ; CATEGORY:
  324. ;    Statistics.
  325. ;
  326. ; CALLING SEQUENCE: 
  327. ;    KMEANS, Data, Cluster, VarName, CaseName
  328. ;
  329. ; INPUTS: 
  330. ;       Data =  a (C,R) dimensioned array where R is
  331. ;        the number of cases to be partitioned and C is the number 
  332. ;        of variables.
  333. ;
  334. ; KEYWORDS:
  335. ;    Number = the number of clusters in final
  336. ;                  partition. The default is 2
  337. ;    Iter =   the number of iterations used to assign
  338. ;                  cases to clusters. The default is 50.
  339. ;    Norm =   flag, if set, to signal whether to
  340. ;                  normalize the variable values in Data.
  341. ;                  Values normalized only if Norm=1. 
  342. ;
  343. ;    Missing = missing data value. If undefined,
  344. ;                   assume no missing data
  345. ;
  346. ;    List_Name= name of output file. Default is to
  347. ;                    the screen.
  348. ;    ClustMeans=array of cluster means. If defined on
  349. ;                    entry, CM(i,j) = mean of ith variable
  350. ;                     in j th cluster on exit
  351. ;    ClustSTD = array of cluster standard deviations.
  352. ;                    If defined on entry, CS(i,j) =
  353. ;                     standard deviation of ith variable
  354. ;                    in jth cluster on exit
  355. ;    SSBetween = array of sum of squares  between
  356. ;                    clusters. If defined on entry,
  357. ;                    SSB(i)= sum of squares for ith
  358. ;                    variable.
  359. ;    SSWithin = array of sum of squares within clusters.
  360. ;                   SSW is the  analogue for SSB
  361. ;    VarName= one dimensional array of C variable names
  362. ;    CaseName= one dimensional array of R case names
  363. ;
  364. ; OUTPUT:
  365. ;    Summary statistics for each variable for each cluster, overall
  366. ;    analysis of variance for each variable.
  367. ;
  368. ; OPTIONAL OUTPUT PARAMETERS: 
  369. ;    Cluster    a one dimensional array. Cluster(i) = the cluster containing
  370. ;    case i in the final partition
  371. ;
  372. ; RESTRICTIONS:
  373. ;    None.
  374. ;
  375. ; COMMON BLOCKS:
  376. ;    None.
  377. ;
  378. ; SIDE EFFECTS:
  379. ;    None.
  380. ;
  381. ; PROCEDURE:
  382. ;    Adapted from algorithm in Clustering Algorithms by Hartigan, 
  383. ;    Wiley Series in Probablity and Mathematical Statistics, Chapt.4.
  384. ;    Function kmeans1 implements a function that given a partition P
  385. ;    returns a partition P' in the same neighborhood with reduced in group
  386. ;    error.
  387. ;
  388. ;    Function is called repeatedly until it finds a fixed point or local
  389. ;    minimum. Kmeans1 recomputes cluster means after each reassignment.
  390. ;           
  391. ;    Procedure Kmeans successively finds partitions with the starting 
  392. ;    partition for K the final partition for K-1 with the case farthest
  393. ;    from its cluster mean split off to form a new cluster.
  394. ;-
  395.  
  396. ;on_error,2
  397.  Common KBlock,B
  398.  
  399.  SD=Size(Data1)
  400.  
  401.  
  402.  if(N_ELements(LN) EQ 0)THEN  Unit = -1  $
  403.   else openw,unit,/Get,LN
  404.  if(SD(0) NE 2) THEN BEGIN
  405.    printf,unit,  $
  406.    " kmeans- Data must be a two dimensional array."
  407.    goto,Done1
  408.  ENDIF
  409.  
  410.  Data = Data1 
  411.    
  412.  if (N_Elements(M) EQ 1) THEN BEGIN          
  413.        Data=ListWise(Data,M)
  414.        if N_elements(Data) le 1 THEN $
  415.        if(Data EQ -1) THEN BEGIN
  416.          printf,unit,   $
  417.      "kmeans- halted since all cases have missing data."
  418.          goto,DONE1
  419.        ENDIF
  420.        SD = size(Data)
  421.  ENDIF
  422.  
  423.  
  424.   
  425.  C=Sd(1)
  426.  if SD(0) eq 1 then R = 1 else R=SD(2)
  427.  
  428.  
  429.  if( N_Elements(N) EQ 1) then if(n GT R) THEN BEGIN
  430.    printf,unit, $
  431.  "kmeans- Partition number must not be greater than number of cases"
  432.     goto,DONE1
  433.  ENDIF else PN=n else PN= 2
  434.  
  435.  if(N_Elements(IT) NE 0) THEN Iter= IT else Iter =50
  436.  
  437.  if(N_Elements(nr) NE 0) THEN  if ( nr EQ 1) THEN BEGIN
  438.         D=Data
  439.         Data= Normal1(Data,R,C)  
  440.   ENDIF
  441.  
  442.  
  443.  
  444.  
  445.  
  446.   Cluster = Fltarr(R)    ; Clus(i)= cluster containing case i
  447.                          ; initially all cases in same cluster
  448.   Temp=Replicate(1.,C)
  449.  
  450.   for i=0l,PN-1 DO BEGIN   
  451.                   ; successively, construct cluster partitions
  452.       kmeans1,data,cluster,iter,i+1
  453.       max2=0
  454.       ind1 = -1
  455.       for j=0l,i do begin
  456.                     ;Find case farthest from cluster mean 
  457.           Ind= where(cluster EQ j,count)
  458.           if(count NE 0) THEN BEGIN
  459.              A=Data(*,IND)-B(*,j)#Replicate(1.,count)
  460.              if count eq 1 THEN BEGIN
  461.                 here = 0
  462.                 max1 = Total(Temp *A^2)
  463.              ENDIF ELSE BEGIN
  464.                max1= max(Temp # (A*A),here)
  465.                here=Ind(here)
  466.             ENDELSE
  467.              if(max1 GT max2) THEN BEGIN
  468.                ind1=here
  469.                max2=max1
  470.              ENDIF
  471.             
  472.           ENDIF
  473.       ENDFOR
  474.  
  475.     if(i NE PN-1 and ind1 ne -1) THEN Cluster(ind1)= i+1
  476.  
  477.   ENDFOR
  478.  
  479.  if (N_Elements(NP) EQ 0) THEN                          $
  480.    OutPut1,cluster,Data,CaseName,VarName,unit,SD1,SB1,SW1
  481.  
  482.  if(N_Elements(CS) NE 0) THEN CS=SD1   
  483.  if(N_Elements(SB) NE 0) THEN SB=SB1   
  484.  if(N_Elements(SW) NE 0) THEN SW=SW1   
  485.  if(N_Elements(CM) NE 0)  THEN CM=B
  486.  
  487. DONE1:
  488.  if( unit NE -1) Then Free_Lun ,unit    
  489.  RETURN
  490.  END                        ; kmeans
  491.